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Abstract 

A semiclassical model of charge transport in a semiconductor superlattice 
is solved, using moments in the wavenumber direction and finite elements 
in the spatial direction (first order). The selection of numerical methods 
guarantees the conservation of current while allowing for high accuracy re- 
sults. When a dc voltage bias is held between the ends of the sample, 
self-sustaining oscillations of the current through the superlattice are ob- 
served in a narrow range of voltages, the calculated solution displayed the 
expected accuracy: Spectral convergence in the number of moments used, 
and first-order convergence in the number of grid-cells. This result paves 
the way for higher-order methods (in the spatial direction) and the numeri- 
cal solution of more complex models of charge transport including quantum 
models based on the Wigner function. 

Keywords: Semiconductor superlattice, kinetic equation of 
Boltzmann-Poisson type, contact boundary conditions, self-sustained 
current oscillations, spectral methods 
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1. Introduction 

Bloch oscillations are coherent oscillations of the electron position inside 
an energy band of a crystal under an applied electric field. Their frequency 
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is proportional to the field times the lattice constant and it should be larger 
than the inverse scattering time for the oscillations to persist. The neces- 
sary electric field is too large for natural crystals and thus Esaki and Tsu 
suggested in 1970 to construct an artificial crystal with a larger effective lat- 
tice constant called a superlattice (SL) |ET70| . The simplest SL example is 
formed by epitaxially growing many identical periods comprising a number 
of layers of two different semiconductor materials |BG05j . The difference 
in the energy gaps of the component semiconductors causes the conduction 
band of the superlattice to be a periodic succession of barriers and wells with 
typical periods of several nanometers. Provided the lateral extension of a 
SL is much larger than its period, it is a quasi one-dime nsional (ID) sys- 
tem. Damped Bloch oscillations of terahertz frequency were first observed 
in 1992 in such undoped semiconductor SLs whose initial state was pre- 



pared optically FLS + 92 . These SLs had finitely many spatial periods and 
were subject to an appropriate DC voltage bias. In SLs made out of doped 
semiconductors, scattering usually destroys the Bloch oscillations but, in 
theory, they can persist even in the hydrodynamic regime for a SL with long 
scattering times |BAC11| IACBT2] . 

Except for a narrow parameter range, Bloch oscillations are not sta- 
ble states in doped SLs [BAClTJ IACB 12]. However there are other stable 
self-sustained oscillations (SSCO) of the current that are observed in a DC 
voltage biased SL. These oscillations have frequencies in the gigahertz range 
and are caused by repeated formation of electric field pulses at the injecting 
contact of the SL that move forward and disappear at the receiving contact. 
They have been observed in experiments with GaAS/AlAs SL (and with 
other SL based on III-V semiconductors) since 1996 and are the basis of fast 
oscillator devices |HGS + 96| , which have important applications in industry. 

At the most fundamental level, nonlinear transport in SLs has been 
modeled using quantum kinetic equations based on nonequilibrium Green 
functions |Wac02 , Wannier-Stark distribution functions [SRD02 or Wigner- 
Poisson equations {BE051 ABU) . In the latt er case, reduced nonlocal drift- 
diffusion equations for the electric field and the electron density can be 
derived using the Chapman-Enskog perturbation method [BE05, AB10]. In 
the semi classical limit, these equations coincide with those similarly derived 
for semiclassical Boltzmann-type equations [BEP03]. Mathematical models 
at the level of semiclassical kinetic theory go back to the 1970s [KSS72] but, 
in the early work, their analysis was based on simplified reduced rate equa- 
tions (ordinary differential equations) |IS87| IIDS91] which typically ignore 
space-charge effects. Electron transport in a single-miniband SL can be de- 
scribed by a semiclassical kinetic equation couple d to a Poisson equation 
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approximately describing the electric potential due to the other electrons 
[BEPQ3j. The electron density in the x — k space (position and momen- 
tum) satisfies a two-dimensional, non-linear, hyperbolic PDE, coupled to a 
Poisson equation which depends on the average charge density. 

Recently, Cebrian et al jCBC09) numerically solved the kinetic equation 
using a direct approach and showed that self-oscillations are among its so- 
lutions and also studied the relation between these solutions and those of 
the limiting drift-diffusion equation. Their numerical solution was based 
on a hybrid particle/fixed-grid method that used the particles to solve the 
advective terms and used the grid for calculating the solution to the Poisson 
problem and for evaluating the effect of the source term. This has several 
disadvantages: First, it adds a layer of complication as the solution needs 
to be continually projected back and forth from the particles to the grid; 
second, despite the conservative nature of the equations it has not been 
shown that the resulting method is conservative. In fact, due to the use of 
averaging for obtaining point- wise values, there is reason to believe that it 
is not; lastly, it would be difficult to extend this solution method to allow 
for a time-dependent bias voltage and for solving quantum kinetic equations 
away from the semi classical limit. 

Here, we solved the charge-transport equation using moments (Fourier 
basis) in k. This approach has several advantages: 

• The 2-D PDE is transformed into a system of 1-D conservation laws 
(for the coefficients of the moments) which can be solved using a stan- 
dard method (upwind, Godunov method); this is indeed how we solve 
it. 

• The zero moment's (average charge density) equation has no source 
terms, which makes guaranteeing a conservative solution much easier. 

• Since the total density is one of the dependent variables, solving the 
Poisson equation involves only a simple linear equation. 

• The first two moments are the current density and energy density, 
quantities of high physical significance and importance. 

• The resulting method can be generalized to solve more realistic models, 
for example using the Wigner-Poisson quantum kinetic equation |BE05, 
AB10], a more complete collision model [BACH or a time-dependent 
voltage bias. 

• Due to the spectral convergence of Fourier expansion, the computa- 
tional cost is lower using this method, for the same accuracy. 
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Our use of Moment methods is similar to the use thereof in the problem 
of Radiative Transfer. This was first derived formally by Gelbard |Gel60, 
IGel6H IGel62j and has been used by many since then, for example by Frank 
et al. |FLSl"H IFKLY07] . The fundamental idea is to write the solution using 
a family of rapidly converging basis functions and then derive the equations 
for the coefficients. In the case of the Boltzmann equation for rarefied gas 
dynamics, spectral methods have been used after constraining the velocities 
(equivalently, wave numbers) to take values on a bounded domain with 
periodic boundary conditions and modifying accordingly the collision term 
[PROO, FMP06]. In our case, the wave number takes values on a bounded 
interval and the distribution function is periodic in it, so that we do not 
have this additional source of numerical error. 

Inevitably, a numerical approximation will require the truncation of the 
series of basis functions to a finite sum, but the equations describing the 
evolution of the coefficients will not be "closed", that is, it will involve one 
(or more) of the truncated coefficients. An external argument is normally 
needed in order to "close" the resulting equations, and while the zero-closure 
(assume that the truncated coefficients vanish) is easy to implement and 
usually good enough (due to the spectral convergence), other closures, such 
as maximum entropy |Jay57| or optimal prediction |FS11| , can be considered 
as they can lead to significant increase in the accuracy for a given number 
of kept moments. 

The paper is organized as follows. In Section [2] the non-dimensional 
equations that describe the system are presented; in Section [3] the method 
of moments is described as it applies to the equations at hand; in Section |4j 
the implementation of the numerical method is described. Results and con- 
clusions are presented in sections [5] and [6j Nomenclature can be found in 
Section [U 



2. Non-Dimensional Model 

As others before [CBC09, BG05], we non-dimensionalize the charge trans- 
port equations using units which, like all other symbols in this paper, can 
be found in Section HI 

The non-dimensional of equations for the electron density f(k, x, t) is: 

ft + 2^sm(k)f x + ^F(x)f k = - \f FD (k, fx) - f(k)(l + M) + Mf(-k) 
r] r] L 
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Here, / and F are coupled via a Poisson equation for the potential V: 

V xx = F x = n-l, V(0)=0, V(L)=<j)L (2) 

1 - If* 

n{x,t) = -=f Q (x,t) = -= / cos(jk)f(k,x,t)dk. (3) 

where /o is the th Fourier modej^of /. The others modes are given by: 

1 [* 

fj = —= \ cos(jk)f(k,x,t)dk for j >0, and (4) 
/-j = -7= / sin(jk)f{k, x, t) dk for - j < , (5) 

V 71 " ^-7T 

while / FD is the Fermi-Dirac distribution, given by 

f FD (k,n) =alog(l + exp[/i-*(l-cos(A;))]). (6) 
In the definition of / FD , [i = n(n) is the unique value for which 

n=-^/ FD (^) (7) 

Where the Fourier modes of / FD are defined equivalently to those of /. 

2.1. Boundary conditions and initial conditions 

To make the problem well-posed, we need to supply initial conditions for 
f(x, k, 0) and boundary conditions for /(0, k, t) and f(L, k, t). Both require 
one more definition. 

A steady-state solution at a constant field, F, would have vanishing time- 
and space-derivatives. We define this distribution as and use it both in 
the initial conditions and in the boundary conditions. Setting the time- and 
space-derivatives to zero in (jlj, we get a (non-local) ODE for f(°'(k;F,n): 

(1 + M)f^(k) - Mf^(-k) + T e Fd k fM = /FD (fc) /x(n))) (g) 

As shown below, a solution to this equation is straight-forward using Fourier 
series. We assume that the initial condition solution is /(°) with n = 1 and 
F = <p: 

f(x,k,0)=f^(k;ct>,l). (9) 

As for the boundary conditions, we expect that the x = terminal will be 
injecting electrons, and the x = L terminal will be collecting them. We 
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therefore follow others in using a "top-down" approach and require that the 
current at the injecting terminal obeys (dimensional) Ohm's Law j = aF, 
while at the collecting terminal we simply require that the (dimensional) 
electron density is Np. The non-dimensional versions of these BC are: 

j(0,t) = 2fcF{0), n(L,t) = l. (10) 

Here, j is the local current density: 

j(x,t)=q sin(k)f(k,x,t)dk = y/Trsf-i(x,t) (11) 



There are several ways to achieve these requirements. Since Eq. ([I]) is hy- 
perbolic, we may only set boundary conditions where the characteristics are 
going into the domain, that is, for x = we should only set conditions for 
k > 0, and for x = L, we should only set conditions for k < 0. Since we have 
only one condition for every boundary, the problem is under-determined. We 
deviate slightly from the choice made in [CBC09 and use a multiple of /(°) 
as the boundary condition: 

/( °' " > °< f) = jyainW/Swdfc ( 2PF ~ L S ^ k)m *' l) dk ) (12) 
f(L, k<0,t)= r0 f ^} k \ ~ - r f(L, k, t) dk) (13) 



f° T f(°)(k)dk\ Jo 

A quick check shows that with these definitions the BC at x = and x = L 
are satisfied. In the k— direction we impose periodic boundary conditions. 

In summary, the problem consists of advection PDE ([!]) coupled with 
the Poisson problem initial condition and boundary conditions (12 



13) in the x— direction, and periodicity in the k— direction. 



3. Method of moments 



We solve this model using Fourier modes in the A;— direction and Gudunov 
method with wave-splitting on a regular grid in the x— direction. This has 
the advantage of imitating the charge density conservation property that the 
original equations have, and provides us with the important physical vari- 
ables (charge density, current density and energy density) for "free", without 
the need to calculate them from the solution. 

Multiplying Eq. ([T]) by ^sin(jA;), -^cos(jk), or ^= and integrating 
from — 7r to 7T (with respect to k) results in the following system of equations 
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for fj, the Fourier coefficients of /: 

dtf) + ^<)Af , . - / + JTj ^ r) f , = l - (// D - ./•) for j > 2 

dth + ^dj^ = 1 (// D -/ i ) farj = l 

djo + iKyfidJ-x = 1 f / FD - /o) = for j = 



dtf-i + nsd x (V2f -f 2 ) = _(i + 2M)^ forj = -l 

r\ rj 

dtfj + iKdvif-j-i - f-j+i) - J -^^f-j =-(l + 2M)^ for j < -2 

We have used integration by parts, and the symmetry with respect to k of 
the Fermi-Dirac distribution, / FD . 

One advantage of using Fourier moments is that the equation for /o has 
no source term. This encodes the fact that electrons are not created or 
destroyed, they are only moved from one place to another with a current 
density By solving the system with a conservative numerical method, 
we are guaranteed that charge is conserved. 

These equations are exact so long as we take all of the infinite moments 
involved. But, of course, when implementing this method we can keep only a 
finite number of moments, therefore in the extremal equations — for and 
/_7V — there will be missing terms: f~N~i for Jn and /jv+i for f~N- This 
is a standard problem in moment methods, the solution thereof is referred 
to as "moment closure". In this paper we use what is known as the "P/v 
closure" which assumes that the missing moments vanish. Other closures 
may be pursued at a later time. 

For a finite number of moments, the moment equations form an advection- 
reaction PDE (in t and x) for the vector of moments f : 

% + irUt = ^{T e F(x)S 1 i + 5 2 f + f FD (/i(n(^)))}, (14) 

where A is the advection matrix given by (15), f FD is the vector of Fourier 
coefficients of / FD , and Si &: S2 are matrices given by (16). 

If we take N positive and N negative moments (and a zero moment), we 
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can write the of (14) as follows: 



( fN \ 

h 

A 
/o 
f-i 



f-N+l 
V f-N ) t 



+ 7T? 
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1 

h 

h 

b 
f-i 







(15) 



The dashed lines separate the negative moments from the non-negative ones 
as a visual aid. With the same notation for the vector of values, the 
matrices Si, S2 are given by 



/ 



Si 



1 



\~ N ■ / V i 

As mentioned earlier, using these moment the Fourier transform of /(°) 



N 



( 



So 



1 



1 + 2M 



1 + 2M j 
(16) 



from Eq. ([8j) is easy to find, as it must solve Eq. ( 14 ) with both derivative 
terms omitted: 

= TeFSxfW + S 2 f(°) + f FD (/x(n)). (17) 



In other words, 



f(0)( n;F ) = -( Te FS 1 + S 2 )- 1 f FD (/x(n)). (18) 



Thus the initial condition is easy to write in Fourier, what about the bound- 



ary conditions? Since the boundary conditions (12) and (13) depend on 



Fourier modes of the truncated density function 0(— k)f(0, k) and Q(k)f(L, k), 
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we reconstruct / at the boundaries from the Fourier coefficients, truncate 
the resulting function as appropriate, and calculate the Fourier coefficients 
of the truncated function. While this isn't very efficient, it only has to be 
done at the boundaries, and therefore relatively cheap, computationally. 

4. Implementation 

The complete problem consisting of linear advection, source terms, and 
the coupled Poisson equation are solved using operator splitting, alternating 
between an advection step and a source step. The electric field is calculated 
from the solution of the Poisson equation before it is needed in the source 
term and boundary conditions. The source consists of three terms that can 
be more accurately solved separately than together. Therefore, we also use 
operator-splitting for the source step itself 

4-1. Advection 

To solve the linear advection system, we use wave-splitting following 
LeVeque's book|LcV02 . This means that the numerical values represent 
cell averages and at every time-step, we write the difference between neigh- 
boring cells as a sum of eigenvectors of the advection matrix, and calculate 
the change to the cell-averages due to upwind advection of these wavea^j 
This is a first-order approximation that is consistent with the conservation 
properties of the problem, and is therefore guaranteed to conserve the charge 
density. 

Given a Riemann problem initial condition (constant solution at each 
cell), each eigenvector of the matrix A corresponds to a "wave" that travels 
at a specific speed A (the corresponding eigenvalue). These waves can be 
followed as they travel forward (A > 0) or backwards (A < 0) and the cell 
averages can be adjusted accordingly: 

£m+l £ m 1 At 



h 



>(fr-f^i)+^-(f^i-fr)] (19) 

Where A^ are the right- and left-going advection velocities given by 

A + = R(A) + R-\ A" = RiAyR- 1 . (20) 



The time-step is chosen small enough so that waves from neighboring cells cannot 
interact 
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For stability we keep ^ smaller than the fastest wave in the system. As 
iV— >oo, the largest eigenvalue of A approaches 1 (this encodes the maximal 
value of sin(fc)) thus, CFL condition we use 



At = .95 



h 

7IX 



(21) 



At the boundaries, we need to provide a value of the solution outside 
the domain. Since only the waves that go into the domain affect it, we can 
provide the values on both the ingoing and outgoing parts of the solution 
and let the up-winding take care of moving the information in the correct 
direction. 



Following Eq. (12), for the left boundary, x = 0, we let a "ghost" cell 
have the value 



~ _ f(°) (2/3F 



(22) 



Similarly, at the right boundary, x = L, we define another "ghost" cell with 
the solution 



f(o) 



I N X +1 



T[Q(-k)fW(k)] 
The inverse Fourier operator is defined as: 



(23) 



g(k) = f-'Mk) = + ^ r J29-jMjk) + 9 j cos(jk). (24) 
4-2. Sources 

There are three sources terms to contend with and our ability to solve 
them analytically differs between them. The linear ones are each trivial 
to solve exactly as the eigenvalues of Si and S2 can be calculated once in 
advance. Using a spectral decomposition of a matrix S: 

RDR- 1 = S, where D is a diagonal matrix of eigenvalues (25) 

The solution to 



1 



N 



f t = Si, with f(t) given 



(26) 



can be written as 



f (t + At) = Rexp(AtD)R- L f(t). 



(27) 
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The source term that comes from the Fermi-Dirac distribution is non- 
linear but while it only depends on /o, the resulting term is non-zero only 
for positive j terms. This implies that it can be solved with a first-order 
integrator — we use Forward Euler — with no loss of accuracy. 

The three sources are put together using first-order Operator-Splitting: 
taking first a step with the non-linear term, then with S± and finally with 
$2- Since the advection is computed to first-order, any effort for calculating 
a second-order solution of the source would be mostly lost. 

4-3. Fermi-Dirac Distribution 

One of the main bottlenecks of the previous papers was in the inversion 
of the chemical-potential function fj,(n). The problem is that this inversion 
is needed at every time-step, at every grid-point X{. A relatively fast solver 
can be written using Newton's method, but even so, the sheer amount of 
calculating is time-consuming. We have elected to use pre-calculation and 
then evaluation using a piecewise cubic Hermite interpolant (for fi) and a 
spline (for the Fourier coefficients of / FD ). By calculating once the value of 
\x (and also f FD ) at sufficient values of n between and Erl we can guarantee 
that the error of the interpolation is smaller than a desired accuracy. The 
evaluation time is a fraction of that when using Newton's method. 

The error of the Hermite interpolation of [i{n) where n S [n^nj+i] is 
bounded by 

Error < II^V)IUKn t+l] (n _ nj)2(n _ nj+i)2 (2g) 

For small values of \x the function \i{n) behaves asymptotically like log(n) 
and this is also where the large fourth derivatives are found. We therefore 
estimate the fourth derivative of [J,(n) for re < 1 as 

/u (4) (n) ss for n < 1 (29) 

We use this to estimate the accuracy of our approximation of \i and when it 
is not accurate enough, we perform a few Newton iterations until the desired 
accuracy is achieved (starting, of course, with the interpolation result). Once 
the correct values are found, the Hermit interplant is modified by dividing 
the offending segment into two parts, thus greatly increasing the accuracy 



3 This is a conservative estimate of the possible range of values we will encounter in the 
simulation. 
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of the interpolation on it. This process was continued until the error was 
less than 10~ 12 for the function /J.(n). 

To evaluate / FD (/i) we use a spline interpolant on 10000 points between 
and 10. 



4-4- Poisson equation 

The solution to the Poisson equation is required in order to find the 
electric field, F, needed in equation Q and in boundary conditions (12 13). 



When using the moment methods, the field will be required in the center of 
each field, to coincide with the other variables. As the field is a first deriva- 
tive of the electric potential (the solution of the Poisson equation) it would 
be optimal if the potential were given on the boundaries of the cells rather 
than the centers. This would also integrate the boundary conditions of the 
Poisson problem most easily, since they too are given on the boundaries. 
The only problem with this approach is that we have N x equations (one 
for each cell) and only N x — 1 values with which to satisfy them (from the 
values of the potential at the internal edges). To find a good candidate for 
this over-constrained problem, we use a finite element approach. 

Our problem therefore is to find the electric potential, V(x), and electric 
field, F, by solving Q: 

V xx = n-1, V(0)=0, V{L)=4)L (30) 

for a given n(x) and (p. Since our first-order advection solver assumes a 
piecewise constant solution space, n(x) is assumed constant within each 
x-cell. We can therefore write it as 



N x 



< hx 



n(x) =]>>,!;, Ux) = \ 1 * " w V" - "~ (31) 
i=1 I otherwise, 

where h is the resolution of the x-grid: h = L/N x , and m is the constant 
value of the charge density in the i— th cell. 

We separate the solution into two parts, one that satisfies the boundary 
conditions and the homogeneous equation V xx = and another that satisfies 
the inhomogeneous equation, but homogeneous BC. The first is, of course, 
x<p/L. It is the second part for which finite elements are used. The elements 
we use for the potential are triangular tpj with j = l...N x — 1: 



if \x — jh\ > hx 

ih\ 

k x h 

To find the field F, we use the following program: 



^■(*) = 1 \*- jh \ ' (32) 
1 1 — - — — otherwise. 
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1. Write the potential V as linear combination of the vy. 



V(x) 



N x -l 

1=1 



(33) 



This gives us N x — 1 coefficients, m, such that for any choice the resulting 
potential satisfies the boundary conditions V(0) = V(L) = 0. 



2. Write the Poisson equation (30) using the density as in (31) and the 
potential as in (33). Multiply by <pk and integrate by parts: 



[L*T*-1 . L /n x 

/ Y] -Vitpi, x (x)ipk,x(x) dx = / y7n;-l)l 



(x) )ijj k) x(x)dx. 



(34) 



3. Switch the order of the integral and the sum on both sides of the inequal- 
ity. This results in a simple linear equation whose solution is the finite 
elements solution to the Poisson problem: 



Sv = M(n - 1) 



(35) 



The vectors v and n are the vectors of coefficients Vi and Uj, and the two 
matrices, S and M are given by: 



S 



( 



V 



■l 



-1 2 



/ 



M 



( 



\ 



1 1 



1 1 



/ 



(36) 



Here S is a square matrix of size (N x — 1), while M has N x — 1 rows, 
and N x columns (implemented as sparse matrices). These matrices arise 
from the integrals in ( |34[ ) The complete solution (including the boundary 
conditions) can thus be written as 



where Xi 



v = x(j)/L- S~ 1 M{n-i). 
ih and 1 is a vector of ones. 



(37) 



4. To calculate the field, F, we numerically differentiate the potential, V, 
thus obtaining the field in the middle of each cell. This is what we need 
for the source term. At the boundaries, we extrapolate from the two 
closest values of F. 
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The result is equivalent to solving using a simple divided differences ap- 
proach, with the effective density n at each edge equal to the average of the 
two densities in the neighboring cells. 



5. Results 

The method described in the previous section exhibits self-sustained os- 
cillations very similar to the ones found in the Cebrian paper |CBC09] . for 
average bias field 4> > 1 (see Fig. [2j. Specifically, we see that the result- 
ing charge distributions are similar, and that the resulting current densities 
shows self-sustained oscillations that have similar temporal frequency and 
similar range. 



We calculate the error by comparing j, the current density (11), ob- 
tained at different choices of N x and N m with one obtained at N m = 15 
and N x = 177827. The I 2 distance between the current densities is the re- 
ported error. Both subfigures in Fig. [6] show that the method displays the 
expected convergence: First-order convergence in the number of x-cells and 
spectral- convergence in the number of moments. However, due to the spec- 
tral convergence as the error due to the truncation of moments is quickly 
over-shadowed by the error proportional to h (as evidenced by the horizontal 
plateau. This implies that a second-order implementation of the advection 
and operator splitting would likely result in higher-order solution in the 
number of x-cells, and therefore with much more accurate results for the 
same computation effort. 

The same figure also evidences an unexpected dependence on N m . It 
is not monotone: the accuracy of N m that are 1 (mod 4) is much lower 
than that of those that are 3 (mod 4) , though both subsequences seem to 
have the same limiting plateau as iV m — t-oo. Other problems that have been 
solved using moment methods also exhibit such non-monotone convergence, 
see for example Davison's book |Dav57j . In those cases the cause of the poor 
accuracy is the existence of a zero eigenvalue in the advection matrix. To 
avoid the lower accuracy approximations other authors used even number 
of moments in their calculations. In this paper, we used odd number of 
moments due to the inherent symmetry of the problem and the physical 
interpretation of the zeroth moment (charge density). Therefore, all our 
simulations have a zero eigenvalue. The difference between 1 (mod 4) and 
3 (mod 4) is not in the existence of the zero eigenvalue; an examination 
of the eigenvectors that correspond to the vanishing eigenvalue shows that 
in the 1 (mod 4) the appropriate eigenvector has a non-vanishing /o term, 
while in the 3 (mod 4) the /o term vanishes. We conjecture that the error 
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Figure 1: The mean current J(t) during the transition and two self-sustaining cycles. The 
period of oscillation is 100£o- The markers indicate locations of "snapshots" shown in the 
following figures. 

accumulation in the /o term due to the zero eigenvector is responsible for 
the higher error in the 1 (mod 4) cases. The fo term is more important than 
the rest since many parts of the problem depend on it. 

We compared the results when using the original boundary conditions 
set out in |CBC09| and found very little resulting difference. This affirms 
the claim that as long as the physical constraints are satisfied (Ohm's law 
at the emitting terminal and charge neutrality at the collecting terminal) 
the specific details of the BC are not very important. 

To study the stability of the system and its response with other values of 
the voltage, (ft, we did two slow (non-dimensional time t = 20000) "sweep" 
with (j) varying from to 4 and back. The resulting mean current is shown 
in Fig. [7] The results show that in the region .92 < cfi < 1.12 the system can 
sustain either a constant solution or a self oscillating one. This behavior 
is consistent with a scenario in which the stable self-oscillations appear as 
a subcritical bifurcation from the stationary state at a critical bias in the 
previous region. For the related drift- diffusion model of the Gunn effect, 
such a scenario is realized when the nondimensional length is large enough 
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Figure 3: The charge density n(x,t) during one period of the solution. 



Electric Field, F 




10 20 30 

Position in superlattice, x 



Figure 4: The electric field F(x, t) during one period of the solution. 
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Figure 5: The energy density during one period of the solution. 
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Figure 6: The accuracy of the method as a function of N x and N m . On the left N m = 3 
(mod 4) and on the right N m = 1 (mod 4). The color indicates log(N x ) as per the color-bar 
on the right (N x was taken to be the integer part of 10 fc ^ 4 with k = 8 . . . 20. 
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Figure 7: The mean current during a voltage sweep. In blue the voltage was increasing, 
and in red — decreasing. 

|BH95, KHB96]. This could have interesting applications, as one may be 
able to encourage the system to pick one behavior over the other by external 
stimulus. 

For (f) < .92 the system does not sustain self-oscillations, while for <p > 
1.12 it not only sustains them, but the constant solution becomes unstable 
(sub-critically). For even larger values of eft (0 ~ 3) the behavior seems 
erratic, and the high derivatives encountered at the collecting terminal put 
in question the validity of the results. We intend to repeat these test with 
a second-order solver to verify. It seems that at these high values of (j) once 
again only the constant solution is possible, but the accuracy of the solution 
deteriorates at such large values of 4> due to the resulting high x— derivatives 
of the solution. The results do not provide a clear determination of whether 
the transition is sub- or super-critical. 

6. Outlook and Conclusions 

This paper is a proof of concept — showing a solution to a model of charge 
transport in a super lattice using moments. We used the simplest approaches 
whenever possible, for example we implemented the integrator with only 
first order accuracy in the x— direction. This allows one to seek improved 
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accuracy in a future study, while knowing that the results, self-sustained 
current oscillations, do not depend on the higher-order method. 

Following are several ways one could improve the accuracy of solution 
and provide a more complete solution to the charge transport problem: 

• Use a second-order solver for the advection, Poisson problem, and 
operator splitting. 

• Use a different moment-closure model. For example one could use the 
high moments of or / FD to close the moment equations. 

• It may be possible to calculate the maximum entropy moment clo- 
sure. This implies finding the most likely distribution given the lower 
moments, and using the moments of that distribution to close the 
equations. 

We have shown that moment methods can be used to solve the problem 
of change transport in a superlattice. The main difficulties of the orig- 
inal problem (namely, the non-local character of the collision kernel and 
the integro-differential character of the Poisson problem) are neatly diffused 
by using a Fourier basis in the k— direction. It provide the charge den- 
sity (required for the Poisson problem) as a dependent variable, enforces 
the periodic boundary condition in k naturally, and cleanly transforms the 
non-local collision term into a simple matrix multiplication. The resulting 
method can be adapted to accommodate other terms that may appear in 
a less primitive model, and could also be improved by smarter integration 
methods and moment closure. 
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8. Nomenclature 



Physical Constants 

eo The permittivity of vacuum, 8.8541 x 10™ 12 C/(m V) 

e The charge of an electron, 1.60218 x 1CP 19 C 

h 1.05457 x 10~ 34 J s. Dirac's constant 

mo The rest mass of an electron, 9.10938 x 10 _31 kg 

u e 9 x 10 12 Hz, the collision frequency of electrons with each other 

v\ 18 x 10 12 Hz, the collision frequency of electrons with impurities 

Derived Constants 

I dw + ds, the period of the superlattice, 4.75nm 

m* Electron's effective mass (0.0674,/ + 0.15d B )^ = 7.64191 x 10~ 32 kg 

Functions 

f(°'(k) short-hand for f(°\k; F, n) with F and n from the context 
f(ty(k; F,n) The steady-state distribution for a given F and n, see eq. (8) 
/J The jth moment of the Fermi-Dirac distribution / FD (fi) 
fj The jth moment of a (non-dimensional) electron density /, see Eqs. 

/ FD (/c;/u) The Fermi-Dirac distribution, see eq. (5) 

fw The Fourier modes of the equilibrium distribution f(°\ see eq. (18) 
P[ • ] The Fourier operator, resulting in a 2N + 1— vector, see eqs. 
F~ 1 [ ■ ] The inverse Fourier operator, resulting in a function of k, see eq. (24) 
Q(k) The Heaviside function Q(k) = 1 for k > and zero otherwise 
Numbers 

a A pre-factor in the definition of / FD . = 0.925115 

(3 The non-dimensional contact conductivity 2 ^^ CT = 0.440331 

5 The non-dimensional energy barrier height, = 29.8402 

At The timestep in the numerical method. 

r] The scaled collision frequency, j^- = 0.476181 

h The size of each grid-cell, h = L/N x 
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L Dimensionless length of superlattice device, N p I/xq=A5 

M The ratio between the collision constants: M = 

N The index j of the fastest oscillating moments in the simulation 

N m The number of Fourier modes used in the numerics N m = 2N + 1 

N x The number of grid-cells in the numerical simulation 

4> The average electric field, unless stated otherwise, we show results 
for <fr = 1 

? Non-dimensional advection coefficient. 47r ^ = 0.582189 

M The value such that / FD (M) = 1, 7.10491 
Problem Parameters 

A The difference in base energy between the two materials, 72meV. 

ds The width of the "barrier" material, 3.64nm 

dw The width of the "well" material, 0.93nm 

e r The relative permittivity of the semi-conductor material, 12.85 

Nr, The density of impurities in the semiconductors, 4.57 x 10 14 m~ 2 

Np The number of periods in the superlattice, 157 

a The contact conductivity 250(m fJ) -1 
Scaling Units 

F M The scaling unit of the electric field V^K+g = 2.24519 xl0 6 J/(C m) 

jo The scaling unit of current density, eVM { ND = 1.094761 x 10 9 A/m 2 

t The scaling unit of time, ^ = 0.233338ps 

vm Scaling unit of electron drift velocity, — — - = 68.3296km/s 

xo The scaling unit of x in the superlattice, ^fj^T = 15.9439nm 
Super- and Sub-Scripts 

( • ) _ min(0, • ) 
( • )+ max(0, • ) 

i Subscript denoting the cell in the numerical method, i = 1 cor- 
respond to the first cell, with left boundary at x = 0. i = N x 
corresponds to the last cell, with right boundary at x = L 
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j Subscript denoting the Fourier mode, see eq. (4) 

m Subscript denoting the timestep in the numerical method 

Variables 

A The advection matrix of the Fourier moments, see eq. (14) 

F The electric field at a point x in the sample, see eq. (2) 

/ The distribution of electrons as a function of x, k, and t, see eq. (1) 

f The vector of Fourier terms of / 

f™ The Fourier terms of / at the i— th cell, during the nth time-step 

J The mean current density: J = x Jo j( x -> i)dx = Jq y/ir(pf-i(x, t) dx 

k The momentum of electrons in the superlattice 

A A diagonal matrix of the eigenvalues of A so that AR = AR 

fi The chemical potential at a given point x, see eq. (6) 

n The total electron density at a point x, see eq. (2) 

R The matrix of eigenvectors of A so that AR = AR 

Si One of the matrices used to generate the source term, see eq. (16) 

S*2 One of the matrices used to generate the source term, see eq. (16) 

t Time, non-dimensionalized with unit to 

V The electric potential in the sample, see eq. (2) 

x The physical dimension of the SL, non-dimensionalized with units xo 
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